Systematic search for low-enthalpy sp 3 carbon allotropes using evolutionary metadynamics 
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We present a systematic search for low-energy metastable superhard carbon allotropes by using the recently 
developed evolutionary metadynamics technique. It is known that cold compression of graphite produces an 
allotrope at 15-20 GPa. Here we look for all low-enthalpy structures accessible from graphite. Starting from 
2H- or 3R-graphite and applying a pressure of 20 GPa, a large variety of intermediate sp s carbon allotropes 
were observed in evolutionary metadynamics simulation. Our calculation not only found all the previous pro- 
posed candidates for 'superhard graphite', but also predicted two allotropes (X-carbon and Y-carbon) showing 
unusual types of 5+7 and 4+8 topologies. These superhard carbon allotropes can be classified into five families 
based on 6 (diamond/lonsdaleite), 5+7 (M- and VK-carbon), 5+7 (X-carbon), 4+8 (bct-C4>, and 4+8 (Y-carbon) 
topologies. This study shows evolutionary metadynamics is a powerful approach both to find the global minima 
and systematically search for low-energy metastable phases reachable from given starting materials. 
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INTRODUCTION 

Carbon can adopt a wide range of structures, from super- 
hard/superdense insulating (diamond, lonsdaleite, hypotheti- 
cal phases hP3, till and tP\2 [1]) to ultrasoft semi-metallic 
(graphite, fullerenes) and even superconducting (doped dia- 
mond [2] and alkali-doped fullerenes [3]). The quest for car- 
bon materials with desired properties is of great interest in 
both fundamental science and advanced technology. One im- 
portant direction in carbon research is the discovery of carbon 
allotropes with advanced mechanical and electronic proper- 
ties. 

It is well known that graphite transforms to the thermody- 
namically stable diamond at high pressures (> 12 GPa) and 
high temperatures (1900-2500 K) |4|. On the contrary, several 
experiments reported that cold compression of graphite pro- 
duces a metastable superhard and transparent phase, clearly 
different from diamond or lonsdaleite, but the exact crystal 
structure could not be determined [5-9]. The difficulty to ex- 
perimentally resolve the crystal structure has stimulated theo- 
retical efforts [10-20 1 . Several structural models were found 
using different techniques. The physical properties of these 
models (M-carbon |TT0][TT]j, W-carbon Q2), 0CI6 (also called 
Z-carbon) 1(1314151. R/P carbon Q3), bct-C 4 EMU) have 
been intensely studied. Simulated x-ray diffraction patterns 
and band gaps of these models are mostly in good agreement 
with experimental data, making it even harder to decide which 
one is the metastable product observed in experiments. On the 
other hand, it is not guaranteed that there is not even a better 
solution for this experimental puzzle. Furthermore, it is likely 
that different metastable phases will be obtained by room- 
temperature compression of different polytypes of graphite, 
or under various non-hydrostatic conditions. This motivates 
us to do a systematic search for low-energy metastable carbon 
allotropes. 

So far, there are several methods to find the ground state 
structures of unknown materials. However, none of them 



are designed to search for metastable states. Our recently 
proposed evolutionary metadynamics method [21] can focus 
on that task. Starting from a reasonable initial crystal struc- 
ture, with this technique one can produce efficiently both the 
ground state and metastable states accessible from that initial 
structure. In this paper, we applied this technique to system- 
atically search for metastable carbon allotropes accessiable 
from graphite. Starting the calculation at 20 GPa from two 
polytypes of graphite (2H and 3R), we easily found the di- 
amond structure (ground state) and a number of low-energy 
metastable structures with sp 3 hybridization which could pos- 
sibly explain 'superhard graphite'. 



METHODOLOGY 

The idea of metadynamics is to introduce a history- 
dependent potential term, which fills the minima in the free 
energy surface so that the system could cross the energy bar- 
riers and undergo phase transitions [22]. This technique is 
usually applied as an extension of molecular dynamics (MD) 
simulation technique [23]. Since it relies on MD simulation to 
equilibrate the system, it often leads to trapping in metastable 
states and amorphization rather than a transition to a stable 
crystal structure. We recently proposed a hybrid method, ba- 
sically a metadynamics-like method driven not by local MD 
sampling, but by efficient global optimization moves lfT0ll24l . 

In this approach [21 ], we start from one known initial struc- 
ture at a given external pressure P. Following Martonak et ah, 
we used the cell vectors matrix hij (also representable as a 
6-dimensional vector K) as a collective variable to distinguish 
the change of state of the system 1 23 1 . For a given system with 
volume V under external pressure P, the derivative of the free 
energy G with respect to h is 



dG 

dha 



= V[h-\P-p)\ 



j» 



(1) 



At each generation (or metastep), many structures are pro- 
duced and relaxed at fixed h, and we select the lowest energy 
structure and compute it internal tensor p. The technique used 
here generates many structures at each metastep, while in tra- 
ditional metadynamics |23|, only one structure is produced at 
each metastep. The cell shape is then updated with a stepping 
parameter 5h 



h im (t+ 1) = h im (t) + 
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where S is the elastic compliance tensor corresponding to 
an elastically isotropic medium with Poisson ratio 0.26, which 
corresponds to the border between brittle and ductile materials 
[ 25 1 and is a good average value to describe both metals and 
insulators. The driving force / = — ^ in Eq. &2t comes from 
a history-dependent Gibbs potential G t where a Gaussian has 
been added to G(h) at every point h(t') already visited in or- 
der to discourage it from being visited again, 



G(t) = G h + J2 We 
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where W is the Gaussian height. Then we compute the vi- 
brational modes for the selected structure according to the dy- 
namical matrix constructed from bond hardness coefficients, 
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Here coefficients a, /3 denote coordinates (x, y, z); coeffi- 
cients a, b, i,j describe the atom in the unit cell; coefficients 
I, m, n describe the unit cell number; ry" is the distance be- 
tween atom i in the unit cell I and atom j in the unit cell n, 
while r 0i '" is corresponding bond distance, and H^ is bond 
hardness coefficient computed from bond distances, covalent 
radii and electronegativities of the atoms [24 1. Note that the 
dynamical matrix corresponds to zero wavevector (extension 
to non-zero wavevectors is straightforward) and unit masses. 
The simulated vibrational modes are used to produce the 
next generation (typically 20-40 softmutated structures). To 
perform softmutation ET1I241 . we move the atoms along the 
eigenvector of the softest calculated mode. One structure 
can be softmutated many times using different non-degenerate 
modes and displacements. The magnitude of the displacement 
(d max ) along the mode eigenvector is an input parameter: 
With relatively small d max and displacements represented by 
a random linear mixture of all mode eigenvectors, the method 
becomes similar to MD-metadynamics in crossing energy bar- 
riers and equilibrating the system. With large d max along the 
softest mode eigenvectors, we obtain the softmutation opera- 
tor 1 24 1, capable of efficiently finding the global energy mini- 
mum. 

The next generation of softmutated structures are produced 
and relaxed in the updated cell. Repeated for a number of gen- 
erations, this computational scheme leads to a series of struc- 
tural transitions and is stopped when the maximum number of 
generations is reached. 



In this work, structure relaxations were done using den- 
sity functional theory (DFT) within the generalized gradi- 
ent approximation (GGA) ll26l in the framework of the all- 
electron projector augmented wave (PAW) [27 ] method as 
implemented in the VASP 11281 |29l code. We used a plane 
wave kinetic energy cutoff of 550 eV for the plane-wave ba- 
sis set and a Brillouin zone sampling resolution of 2ir x 0.08 
A -1 , which showed excellent convergence of the energy dif- 
ferences, stress tensors and structural parameters. 



RESULTS AND DISCUSSIONS 

In a compression experiment at low temperatures (low 
enough to preclude transition to the stable state), the product 
depends on the nature of the starting materials, and on the en- 
ergy landscape (in particular, energy barriers). To fully inves- 
tigate the possible candidate materials, we performed simula- 
tions starting from two different graphite polytypes (graphite- 
2H and 3R), which differ in the stacking of graphene layers. 



Starting from graphite-2H 

We did a preliminary test at 20 GPa starting from the 
graphite-2H structure, and successfully found diamond as the 
ground state, and M-carbon and bct-C4 as metastable states 
lETi In the calculation we set d max =2.5 A, 1^=4000 kbar-A 3 
and (5/i=0.6 A. 

Fig. [TJshows the enthalpy evolution. Graphene layers (Fig. 
Eh), persisted until the 15th generation. Then, upon sufficient 
cell deformation, the layers began to buckle, and the planar 
structure transformed into 3D-networks of ,s/? 3 -hybridized car- 
bon atoms. Lonsdaleite with 6-membered rings (Fig. Eb) ap- 
peared as the best structure in the 16th generation. We ob- 
served in the same generation the bct-C4 structure with 4+8 
membered rings (Fig. El), and M/W-carbon structures con- 
taining 5+7 membered rings (Fig. Ekf), appeared shortly af- 
ter. Lonsdaleite survived for a few generations until it tran- 
formed into a hybride structure made of alternating layers of 
M-carbon and diamond (Fig. Efc, we refer to it as M+D type), 
followed by the transition to another hybride structure made 
of bct-C4 and diamond (Fig. EB, similarly, we refer to it as 
B+D type). Diamond was dominant in the following genera- 
tions. At the 51st generation, the system reverted to graphite. 

The power of the evolutionary metadynamics method lies in 
that it is highly suitable for harvesting low-energy metastable 
structures in addition to the ground state. Those previously 
proposed candidate structures for the product of cold com- 
pression of graphite, bct-C4, M, and W-carbon are all eas- 
ily recovered in a single simulation. More interestingly, we 
also observed many low-energy structures based on 5+7 or 
4+8 topology. The B+D type structure (Fig. EH) observed 
in the simulation is actually the oC 16 structure (sometimes 
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FIG. 1. (Color online) (a) Enthalpy evolution during the compres- 
sion of graphite-2H at 20 GPa; (b) Enthalpy evolution during the 
compression of graphite-3R at 20 GPa (black line: enthalpies for 
best structures with constant cell matrix; red line: enthalpies for best 
structures after full relaxation). 




FIG. 2. (Color online) Structures observed during compression of 
graphite-2H at 20 GPa. (a) graphite-2H (2x2x2 supercell of the 
calculation model); (b) lonsdaleite; (c) bct-C4 with 4+8 membered 
rings; (d) Z-carbon (belonging to B+D type); (e) M-carbon with 5+7 
membered rings; (f) VK-carbon with 5+7 membered rings; (g) M+D 
type carbon. Polygons are highlighted by different colors (quadran- 
gle: turquoise; pentagon: green; hexagon: blue). 



called Z carbon), recently suggested as a candidate for su- 
perhard graphite. Since oC16 inherits layers of bct-C4 and 
diamond, there is no surprise that its thermodynamic proper- 
ties are intermediate between these two structures. The 5+7 
class of structures shows a larger diversity. In some of these 
structures, 5-membered rings form pairs, while in others these 
they are single. The difference of the 5-membered ring pairs' 
orientation leads to two allotropes: M-carbon (Fig. |5£) and 
W-carbon (Fig. [2}). Some structures can be thought of as 
combinations of layers of the M-carbon and diamond struc- 
tures (M+D carbon, as shown in Fig. |2| 



Starting from graphite-3R 



Starting the calculation at 20 GPa from another polytype, 
graphite-3R, which contains three layers per lattice period, 
we again easily found the diamond structure and a number of 
low-energy metastable structures with sp 3 hybridization. Fig. 
plshows the results. Since the initial model has three graphene 
layers, it could form a large variety of M+D and B+D struc- 
tures based on 4+6+8 or 5+6+7 topologies. For instance, we 
observed a B+D structure containing 2 x (4+8) layers and 1x6 
layer (Fig. [3}:), or 1x4+8 layers and 1x6 layer (Fig. [3}l); 
and M+D structure containing 2x(5+7) layers and 1x6 layer 
(Fig. |3^,f). Most strikingly, we also observed another struc- 
ture with a 5+7 topology, which is in Fig. [3^. The projec- 
tions of pentagons and heptagons along the c axis could not 
be separated as in M-carbon, but overlap each other. We ex- 
tracted the 5+7 part from the complex structure, and obtained 
a different configuration with pure 5+7 topology. This crystal 
structure (which we call X-carbon, and the hybrid structure 
from X-carbon and diamond is referred to as X+D type) is 
shown in Fig. HI It is a monoclinic structure with C2/c sym- 
metry, and contains 32 atoms in the conventional cell. We also 
found an unexpected 4+8 topology in an allotrope that we call 
Y-carbon with unique 4+8 ring topology from another sepa- 
rate metadyanmics run. It is an orthorhombic structure with 
Cmca symmetry, containing 16 atoms in the conventional cell. 
The simulated X-ray diffraction patterns of all these structures 
are in good agreement with experimental data (as shown in 
Fig. [5]), suggesting that both X,Y carbon could explain the ex- 
periments on cold compressed graphite. Although all these 
structures show a satisfactory agreement with experimental 
X-ray data, our recent transition path sampling calculations 
lf30ll suggest M-carbon to be kinetically the likeliest product 
of cold compression of graphite-2H. Using other polytypes of 
graphite, or different conditions (non-hydrostatic or dynam- 
ical compression), one might produce alternative allotropes 
found here. Synthesis of these allotropes would be desirable 
in view of their physical properties. 




FIG. 3. (Color online) Structures observed during compression of 
graphite-3R at 20 GPa. (a) graphite-3R (2x2x2 supercell of 
the calculation model); (b) diamond; (c,d) B+D type structures; (e,f) 
M+D type structures; (g) X+D type structures. Polygons are high- 
lighted by different colors (squares: turquoise; pentagons: green; 
hexagons: blue). 
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FIG. 4. (Color online) (a) New allotropes with 5+7 topology, 
X-carbon, space group C2lc, a=5.559 A, b=7.960 A, c=4.752 A, 
f3=l 14.65°. This structure has five non-equivalent Wyckoff posi- 
tions: Cl(0.250, 0.083, 0.949), C2(0.489, 0.809, 0.982), C3(0.000, 
0.200, 0.250), C4(0.247, 0.913, 0.801), C5(0.000, 0.816, 0.250). 
(b) New allotrope with 4+8 topology, /-carbon, space group Cmca, 
a=4.364 A, b=5.057 A, c=4.374 A. This structure has one Wyckoff 
position: C(0.681, 0.635, 0.410). Polygons are highlighted in differ- 
ent colors to show the 5+7 and 4+8 topology. 



Properties 

From evolutionary metadynamics simulations, five fami- 
lies of sp 3 -hybridized structures made by stacking corrugated 
graphene layers and having competitive enthalpies were dis- 
covered: 6 (diamond and lonsdaleite), two classes of 5+7 
topologies (one - M/W-carbon; the other - X-carbon), and two 
classes of 4+8 (one - bct-C4; the other - F-carbon). The en- 
thalpies of different carbon phases as a function of pressure 
are presented in Fig. [6] Apart from the prototypes, we also 
included hybrid structures (lowest enthalpy B+D, M+D and 
X+D carbon, see crystallographic data in Supplementary Ma- 
terials). At elevated pressures, all these allotropes become 
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FIG. 5. Comparison of simulated x-ray diffraction patterns of pro- 
posed models and graphite with experiment [9|. The patterns were 
simulated using Accelrys Materials Studio 4.2 software with the X- 
ray wavelength of 0.3329 A. 



more stable than graphite. For the prototypes, M/W-carbon is 
energetically more favorable than bet- and X-carbon. Lower 
enthalpies are obtained by combining layers of these struc- 
tures with layers of diamond. For the models under consider- 
ation (up to 4 graphene layers), B+D (lx(4+8) + 2x 6 lay- 
ers, Fig. pB) tends to have the lowest enthalpy, while M+D 
(2x(5+7) + 2x6 layers, Fig. |2ji) is quite competitive (only 8 
meV/atom higher than B+D). X+D (Fig. [3ji) is 50 meV/atom 
higher than B+D, indicating that X-carbon has poorest possi- 
bility to interface with diamond. 

We also computed the mechanical and optical properties 
(see Supplementary Materials). Similar to previous theoret- 
ical investigations lHTHT6l[T8l[T9l , all of these candidate al- 
lotropes exhibit high hardnesses BTI and bulk moduli, which 
are comparable with those of diamond. Fig. 17] shows the cal- 
culated total and partial electronic densities of states in both 
systems. It can be clearly seen that 2p states exhibit a larger 
overlap with 2s states in diamond, which makes diamond the 
most stable allotrope among sp 3 forms of carbon. The mag- 
nitude of overlap determines the order of stability: M-carbon 
> X-carbon > bct-C4 > F-carbon. The DFT band gaps of 
M-carbon, X-carbon, bct-C4, and F-carbon are 3.6, 3.8, 2.7 
and 2.9 eV, and we should bear in mind that DFT always un- 
derestimate the band gaps - so the real gaps are larger, and all 
of these sp 3 -allotropes should be transparent colorless insula- 
tors. 



CONCLUSIONS 

In summary, we performed a systematic search for 
metastable allotropes of carbon that can be synthesized by 
cold compression of graphite by using the recently devel- 
oped evolutionary metadynamics technique [21 J. Starting 
from 2H- or 3R-graphite, at 20 GPa we easily found dia- 
mond as the ground state and observed a large variety of 
low-energy metastable sp 3 carbon allotropes accessible from 
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FIG. 6. (Color online) Enthalpies of various carbon structures rela- 
tive to graphite-2H. 
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FIG. 7. (Color online) Total and partial densities of states of different 
carbon allotropes. 



graphite. Apart from diamond, lonsdaleite and their poly- 
types, we summarize four other families of low-enthalpy car- 
bon allotropes which can be obtained by cold compression 
of graphite, (i)5+7 topology (M/W-carbon); (ii)5+7 topology 
(X-carbon); (iii)4+8 topology (bct-C 4 ); (iv)4+8 topology (Y- 
carbon). All of these structures are consistent with experi- 
mental data on 'superhard graphite', and are predicted to have 
excellent mechanical properties, but transition path sampling 
calculations [ 30 1 unequivocally show M-carbon to be the like- 
liest product of cold compression of graphite. Yet, the al- 
lotropes predicted here could be synthesized with a differ- 
ent experimental protocol and starting from different materials 
(graphite-3R, turbostratic graphites, fullerenes, or nanotubes, 
etc). We find it particularly encouraging that all the previously 
proposed structures and two topologically interesting and pre- 
viously not described ones (X- and Y-carbon) were found in 
a systematic way, using just one calculation per starting ma- 
terial (graphite-2H or 3R). Our work shows that evolutionary 
metadynamics is a powerful method for efficiently finding not 
only stable but also low-energy metastable structures. 

Calculations were performed at the supercomputer of Cen- 
ter for Functional Nanomaterials, Brookhaven National Lab- 



* qiang.zhu@stonybrook.edu 
[1] Q. Zhu, A. R. Oganov, M. A. Salva do, P. Pertierra, and A. O. 

Ly akhov, |Phys. Rev. B 83, 193410 (20lT)l 
[2] E. A. Ekimov et al , Nature 428, 542 (2004) 
[3] K. Tanigaki et al, Nature 352, 222 (1991) 

[4] T. Irifune et a/., [Nature 42 1, 599 (2003)[^ 

[5] R. B. Aust and H. G Drickamer, Science 140, 817 (1963) 

[6] M. Hanfland, K. Syassen, and R. Sonnenschein, Phys. Rev. B 

40, 1951 (1989) 

[7] Y. X. Zhao and I. L. Sp ain,|Phys. Rev. B 40, 993 (1989)| 
[8] W. Utsumi and T. Yagi, Science 252, 15 42 (1991)| ~ 

[9] W. L. Mao et al, Science 302, 425 (2003) 

[10] A. R. O ganov arid C. W. Glass, |J. Chem. Phys. 124, 244704| 

(2006)| 

[11] Q. Li et al, [Phys. Rev. Lett. 102, 175506 (2009)| 

[12] J.-T. Wang, C. Chen, and Y. Kawazoe, |Phys. Rev. Lett.^OTl 

075501 (20TT)1 
[13] D. Selli, I. A. Baburin, R. Martonak, and S. Leoni,[Fhys7Rey] 

B 84, 16 1411(2011)1 

[14] Z. Zhao et al, Phys. Rev. Lett. 107, 215502 (2011) 
[15] M. Amsler et al, Phys. Rev. Lett. 108, 065501 (2012) 

[16] H. Niu et al, Phys . Rev. Lett. 108, 135501 (2012)| 

[17] R. Baughman, A. Liu, C. Cui, and P. Schields, Synthetic Metals 

86,2371 (1997)1 
[18] X.-F. Zhou, G.-R. Qian, X. Dong, L. Zhang, Y Tian, and H.-T. 

Wang, Phys. Rev. B 82, 134126 (2010) 
[19] K. Umemoto, R. M. Wentzcovitch, S. Saito, and T. Miyake, 

|Phys. Rev. Lett . 104, 125504 (2010 )] 

[20] V Greshnyakov and E. Belenkov, JETP 113, 86(2011) 

[21] Q. Zhu, A. R. Oganov, and A. O. Lyakhov, CrystEngComm , 

(2012)| 

[22] A. Laio and M. Parrinello, |Proc. Natl. Acad. Sci. 99, 12562| 

(2002)| 
[23] R. Martonak, A. Laio, and M. Parrinello, [Phys. Rev. Lett. 90,| 

075503 (2003) 
[24] A. O. L yakhov, A. R. Oganov, and M. Valle, |Comp. Phys!] 
Comm . 181, 1623 (2010) 

[25] S. F Pugh, Philos. Mag. 45, 823 (1954). 

[26] J. P. Perdew, K. Burke, and M. Ernzerhof, [Phys. Rev. Lett. 77,| 
3865 (1996) 



[27] P. E. Blochl, Phys. Rev. B 50, 17953 (1994) 



[28] G. Kresse and J. Furthmuller, Phys. Rev. B 54, 1 1 169 (1996) 
[29] G. Kresse and D. Joubert, [Phys. Rev. B 5 9, 1758 (1999) 
[30] S. E. Boulfelfel, A. R. Oganov, and S. Leoni, Sci. Rep. (2012), 
accepted. 



[31] A. O. Lyakhov and A. R. Oganov, Phys. Rev. B 84, 092103 
(2011)1 



